Hopf bifurcation in a flux qubit coupled to a nanomechanical oscillator 
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We study the nonlinear semiclassical dynamics of a driven flux qubit coupled to a nanomechanical 
oscillating beam. While both of these systems are dissipative, their mutual coupling can effectively 
suppress the dissipation. This can lead to a loss of stability and to an emergence of synchronized 
self-excited oscillations of the system as a whole, at a time scale set by the mechanical beam. In this 
article we argue this by obtaining a set of semiclassical, nonlinear equations of motion for the two 
coupled subsystems, and showing that this system undergoes a supercritical Hopf bifurcation as the 
coupling is increased. We derive analytical expressions for the critical coupling coefficient and the 
renormalized mechanical dissipation coefficient and frequency, and give a complete characterization 
of the limit cycle behavior when the qubit and the oscillator are near resonance. 

PACS numbers: 85.85.+j, 85.25.Dq, 05.45.-a 



I. INTRODUCTION 

The nano-electromechanical system comprised of a 
driven flux qubit coupled to an oscillating beam has re- 
cently become the subject of a great deal of researcbir— , 
due to its potential for the preparation of macroscopic 

and 
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non-classical states of the nanomechanical bear 
for the realization of high-precision position detectors 
Since the coupling between the qubit and oscillator is 
nonlinear, the semiclassical dynamics of the system may 
exhibit a wide range of typical nonlinear behavior, such 
as loss of stability and limit cycle oscillations. A charac- 
terization of these dynamics is therefore important for an 
experimental realization of these systems, and may also 
hint at possible directions for the preparation of non- 
classical states of the mechanical oscillator. 

In this article we provide a detailed characterization of 
the weakly nonlinear behavior in this system. In partic- 
ular, we show that when the coupling between the qubit 
and the oscillator, which is assumed to be small, is in- 
creased beyond a certain threshold, the system begins to 
exhibit limit cycle behavior, with the qubit and the beam 
undergoing synchronized oscillations at the renormalized 
oscillation frequency of the beam. While the stability 
properties of this system, and of the analogous system 
of a flux qubit coupled to an LC resonator, have been 
investigated in Refs.[5l. [l8l - |23 . our approach is unique by 
the fact that we provide a detailed characterization of 
the dynamical behavior of the system beyond its stabil- 
ity threshold, using analytical tools from the theory of 
nonlinear dynamical systems. Using these methods, we 
derive analytical expressions for the critical coupling co- 
efficients, the renormalized mechanical dissipation coeffi- 
cient, the amplitude and the frequency of the self-excited 
oscillations, in terms of directly measurable experimental 
quantities. We supplement these analytical expressions 
with numerical analysis and simulation. 

The system we investigate consists of a flux-driven rf 
superconducting quantum interference device (SQUID), 
integrated with a mechanical resonator in the shap e of 
a doubly clamped beam. As shown in Refs. |25| - |3C)I . this 



type of system can function as a qubit, exhibiting coher- 
ent quantum dynamics. The motion of the mechanical 
resonator can be approximated as that of a single-mode 
harmonic oscillator. This oscillator affects the dynamics 
of the qubit by varying the amount of flux threading the 
loop of the SQUID, which alters its energy level split- 
ting. The qubit, in turn, affects the dynamics of the 
oscillator by applying a transverse Lorentz force, which 
is proportional to the current in the loop and the mag- 
netic field at the location of the beam. In this article, we 
consider a single Josephson junction (JJ) rf SQUID due 
to its simplicity. In practice, however, the fiux qubit is 
implemented using an asymmetric configuration of three 
or four Josephson junctions known as a persistent current 
qubit^§.i^i^, owing to its superior noise characteristics. 

It has been shown^iiS-^o both theoretically and exper- 
imentally that the coupling between the oscillator and 
the qubit can lead to a renormalization of the oscilla- 
tor's dissipation coefficient through a process known as 
"Sisyphus amplification and cooling" : When the driving 
frequency of the qubit (or one of its harmonics) is larger 
than the energy level separation, the system is said to 
be blue-detuned, and the excess energy in each photon is 
partially transferred to the oscillator, effectively lowering 
its dissipation coeHicient. When the driving frequency is 
smaller than the energy level separation, the system is 
said to be red-detuned, and the opposite occurs^. 

Consider the case of a blue-detuned system. When 
the coupling, which is controlled by the applied magnetic 
field at the location of the beam, is strong enough, the ef- 
fective mechanical dissipation coefficient can become neg- 
ative. When this occurs, the system will undergo a Hopf 
bifurcation and exhibit self-excited oscillations. In this 
article we provide a characterization of these oscillations, 
and qualitatively discuss their consequences for the gen- 
eration of non-classical states of the oscillator. 

The article is organized as follows. In section [ll] we 
derive the semiclassical dynamical equations from first 
principles. This section follows closely the work in Refs. [3 
and 7. In section IIIII we perform an analysis of these dy- 
namical equations. In subsection IIII Bl we determine the 
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stability properties of the linearized system, which are 
identical to those of the full system. In subsection IIIICI 
we derive a complete analytical characterization of the 
limit-cycle behavior of the system when the oscillator 
and the qubit are resonant, which includes expressions 
for the amplitude, frequency and phase of the limit cy- 
cle. Finally, in section ITVl we discuss the possible uses of 
these results and outline directions for future progress. 



II. DERIVATION OF DYNAMICAL SYSTEM 

In this section we derive a set of evolution equations 
for the averages of the operators of the qubit and oscilla- 
tor, under the assumption that these averages factorize. 
We begin with the Hamiltonian of the rf SQUID and 
the mechanical oscillator, from which we derive an effec- 
tive two-level Hamiltonian for the two lowest-lying per- 
sistent current states of the SQUID. We then include the 
effects of the environment by incorporating three reser- 
voirs, which account for the relaxation of the oscillator 
and for the relaxation and dephasing of the qubit. These 
lead to Bloch-Langevin type equations for the qubit and 
a Langevin equation for the oscillator, from which we 
derive the dynamical system under investigation. This 
is found by taking the averages of the equations, and 
assuming that the reservoirs are in thermal equilibrium, 
and that the correlations between the qubit and the oscil- 
lator are sufficiently small to be neglected in the evolution 
equations. 



A. Description of the system 

In this section we recapitulate the main results given 
in Refs. and 0- The system consists of an rf SQUID 
with a single Josephson junction integrated with a doubly 
clamped beam which is free to oscillate-, as shown in 
Fig.[TJ The loop of the SQUID is threaded by a magnetic 
fiux given by 



with 



$t = $^ + L/l + BIX, 



(1) 



where $c is an externally applied flux, L is the self induc- 
tance of the loop, B is the magnetic field at the location 
of the loop, X is the displacement of the beam and / is 
its effective length. 
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Here m is the mass of the beam, C is the capacitance of 
the Josephson junction, Wm is the natural frequency of 
the beam, and = $o^c/27r, with $o = ^/2e and Jc 
is the critical current of the junction. P and Q are the 
momentum of the beam and the charge of the junction 
capacitance, respectively, and we have [X, P] = [$t, Q] = 
ih. 



2. Normalized coordtnates 



Introducing normalized coordinates 
27r$t , 



9 = 
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2ttBI 
27r$e 



-X 
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we write the potential ^ as 



U{(l)t,q)^Ej cos0t~l + — ( 
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+ 2"^ ter 



k- 4>c~ q)' 



(5) 



where /3l = 2ttLIc/^o- In addition, we define a current 
operator 



^ _ 2tt dH 



(6) 



This operator corresponds to the persistent current cir- 
cling the loop. Assuming also that /3l 1, the potential, 
as a function of (pt, has two shallow wells. Each of these 
wells gives rise to a localized phase state, or equivalently, 
to a localized current state. Since the separation of the 
wells is of the order of unity, and Ic is of a macroscopic 
size, these two states correspond to macroscopically dis- 
tinguishable currents, which we label |0) and |0). Fur- 
thermore, the external flux driving is assumed to have 
both DC and AC components, and is given by 

4>c = 4>cO + 0ol COSWdt. 



Hamiltonian 



The combined system of the qubit and oscillator is de- 
scribed by the Hamiltonian Effective Hamiltonian of the qubit and oscillator 

p2 q2 When {(pt) ^ (j)e,q and /3l ^ 1, the two lowest-lying 

energy levels, whose eigenstates are linear combinations 



3 




FIG. 1. (Color online) The qubit-oscillator system, as de- 
scribed in lll Al The inset shows the potential (|3]) as a function 
of (/)t for a certain value of 4>c + g, and both of the localized 
flux wavefunctions with lowest energies. The wavefunction on 
the left-hand side corresponds to the state |0), and the wave- 
function on the right-hand side to lO). The time variation of 
(jic + g in ([3} gives rise to transitions between these two states. 



of |0) and |0), are far from the energy levels above. 
When, in addition, the temperature is sufficiently low, 
it can be assumed that these are the only occupied lev- 
els. In this case, the SQUID can be approximated as a 
qubit, as described in Ref. The circulating current can 
be found using dH]): its eigenvalues are ±/cc, where 



/3l ' 



and (pt.min is the normalized flux at the minima of the 
potential, which can be found approximately by expand- 
ing ([5]) to obtain a quadratic polynomial, giving 



(7) 



The Hamiltonian of the qubit and the oscillator is then 
given by 



So 1 

Hs ^ Y + '^<=i cos WdO + -^l^Ox 




-K/lWin {^0. -f + —hk{a + a^)az 



(8) 



where and cr^ are the Pauli matrices. 




(9) 



and hk = Bllcc y 2mui ^^'^ coupling energy between 
the qubit and the oscillator. 



C. Coupling to the environment 

1. Model of the environment 

The relaxation properties of the qubit and oscillator 
play a pivotal role in the dynamics of this system, and, 
indeed, are the focus of this article. Therefore any valid 
description must include the effects of coupling to the 
environment. We will model the main sources of relax- 
ation as bosonic reservoirs coupled linearly to the qubit 
and oscillator, and consider three types of relaxation pro- 
cesses: thermal relaxation of both the qubit and oscilla- 
tor, and dephasing due to fluctuations in the magnetic 
field threading the loop. The Hamiltonian accounting 
for these reservoirs is 



(10) 



where the boson operators Ci{uj) and cl{uj), with i = 
m, gi, Qip, satisfy the usual commutation relations 

[cj(a;),c](tj')] ^ dijS{uj - uj'), 

where Sij is the Kronecker delta and S{uj^uj') is the delta 
function. 

The interaction between the system and the reservoirs 
is described, in the rotating wave approximation (RWA), 

by 



i^SR 



dwe^'^"a^c™(w) +h.c. (11) 



+h\ ^ I dwe*'^''V+c„i(a;) -Kh.c. 
V 27r 7 

^^V5 / dwe^*'^a,Cg^(w)+h.c., 
where we have used the raising and lowering operators 



1. 

0-+ = -{a^ + i(Ty) 



(12) 



Note that the coupling between the reservoir and system 
operators is assumed to be sufficiently broad-banded and 
slowly varying, on the scale of characteristic frequencies 
of the system, to be a approximated as constant. We have 
also neglected the radiative correction of the oscillation 
frequencies. In addition, we have assumed that the spin- 
flip transitions occur between the two eigenstates of — 
this approximation is valid as long as /lA <^ eg, which is 
the usual case with flux qubits (see, for example, Refs.[28l 
andllil). 
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2. Bloch-Langevin equations 

Solving for the time evolution of the reservoir operators 
in the Heisenberg picture and substituting the result into 
the equations of motion of the system operators, leads to 
the so-called quantum Langevin equation for the oscilla- 
tor and the Bloch-Langevin equations for the qubit. If 
we define envelope operators such that 



cr_(t) =a_(t)e-"'(*), 
a{t) = a(t)e-'"»*, 



(13) 



where 



other system operators, we obtain the following equa- 
tions: 

~ik{a^{t) + a'^^it))a^{t), (18a) 
^f{t) = -71 {(^"^(t) - a,,,q) (18b) 
+ lA (^^(<)e-'''(*) + a5(t)e"'(*)) , 

aR(i) = -|i5«(i) - ^haf{t)e^--\ (18c) 

where 



r]{t) = eo 4>cQt H sinwd^ 



we obtain: 



(14) 



~a.{t) - -ik {ait) + at(t)) a_(<) - (^|- + 7^) ^-(t) 

+ iiAe"'(*V,(i) + t^Vqi{t)<j,{t) (15a) 

aS) = *A (a_(t)e-^''« - a+(t)e*''W) - 7UI + ^40) 
+ 2z {yl{t)d^{t)-d+{t)V,,{t)) , (15b) 

h{t) - -^~a{t) - zifca,(t)e'-'"* + V„(t), (15c) 



with 



71 = 27^2^0 + 1), 
7¥> = 7^(2rio + 1), 

72 = 71 +7v' 

Cz rn = — tanh , 



(19) 



the mean occupation number of the qubit reservoirs. 



D. Passage to rotating frame 



where 
V,i(i) 

V,„(i) 



7]_gi0gl 



47r 



dt^e-*"(*-*°)c,i(to,^), 
dc^e-^'^(*-*°)c,^(to,c^), 
dc^e-^"(*~*°)c„,(to,c^). 



(16) 



The derivation of these equations is standard, and can 
be found, e.g., in complements Civ and Ay of Ref. |H. 
We now average (fTSl) over the degrees of freedom of the 
reservoirs, which we assume are in thermal equilibrium 
with temperature T. Let us denote 



where 



PR 



= TrR [pn(J^] , 

\ k-sT ) 



exp 



TrR, 



exp 



(17) 



is the density operator of the reservoir and fee is Boltz- 
mann's constant, and Tr^ is the partial trace over the 
reservoir coordinates. With a similar definition for the 



We may eliminate the explicit time dependence in 
by performing the RWA as described in appendix C of 
Ref. Is^ . In most experimental setups, Eq is significantly 
larger than fuvd, and the flux driving of the qubit may be 
very strong, witb^i^^ £o4>ci > ^ ^A. We therefore 
need to consider the possibility of multi-photon Rabi os- 
cillations, whereby eo^co — nhujd for a particular integer 
n, and all other harmonics of LOd are far from Eo^gi/^- In 
this case we may define 



= — nujd- 



(20) 



Expansion of e*'''^*^ as a Fourier series then gives, under 
the RWA, 



1——0C 



fUjjA 



(21) 



where J„ is a Bessel function of the first kind. Since all 
terms in (j2ip with I 7^ n will contribute rapidly rotating 
terms to the equations of motion (jl8p , we can perform a 
RWA and neglect them, effectively substituting the factor 



by 



rujJd 
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If we now return to the operators in the non-rotating 
frame, we ehminate the exphcit time dependence, with 
the cost that now the analysis is confined to a muhi- 
photon resonance with a particular n. Note that the 
additional modulation of the energy levels of the qubit 
due to the motion of oscillator is negligible in comparison 
to eo0ci, and has therefore been ignored here. 

E. Semiclassical approximation 

Equations ^TE\\ together with the RWA made in the 
previous section, provide a complete specification of the 
properties of the system, under the assumptions made 
above. If, however, we are only interested in the stability 
properties of the system and in the possibility of existence 
of self-excited oscillations, it is sufficient to consider the 
semiclassical equation formed by averaging these equa- 
tions over the degrees of freedom of the qubit and the 
oscillator. Furthermore, we note that when fc = 0, the 
evolution of the qubit and the oscillator is independent. 
This implies that the covariance function between the 
oscillator and the qubit coordinates is of order fc, which 
allows us to neglect them in the semiclassical equations 
when k is small. Under this assumption, all correlations 
of the system operators factorize, and we arrive at the 
following equations: 

s„(t) = -72S-(i) - iSs_{t) + i^A„s,{t) (22a) 

- ik{a{t) + a*(i))s_(i), 
s.(t) = -7i(s^(i)-a,,cq) (22b) 

+ iA„(s„(t)-s+(t)), 

7m 1 

a{t) = -iojyna{t) ^^i^) ~ "i-^^s^it), (22c) 

where 

and for an operator X, 

{X)s = Trs[pX]. (25) 

Note that we have assumed that the degrees of freedom of 
the reservoir and of the system are statistically indepen- 
dent, and have made the so-called Born approximation 
for the reservoir. 

These equations can also be written in vector form as 

S(i) = A/S(i) + fcg(S(i)) + Seq, (26) 

where M is a 5 x 5 matrix that represents the linear 
part of ([22I1 that does not include the interaction term, 
kg{S(t)) is the interaction term and Soq is the inhomo- 
geneous part of the system. The boldface coordinates 
represent a 5 x 1 vector, with 

S(t) = (s_(t),s+(<),s,(t),a(t),a*(t))T. (27) 



III. ANALYSIS OF DYNAMICAL SYSTEM 

In the case where the system is blue detuned, with 
S = eo4>co/h — nwd < 0, the excess energy of the 
external flux driving is transferred to the mechanical 
oscillatorii2i2£i2^j and is manifested in the system by a re- 
duction of mechanical dissipation. As we will show, when 
the coupling exceeds a threshold, which we denote as kc, 
the mechanical dissipation coefficient becomes negative, 
and the system loses its stability. This loss of stability, 
however, is of the "soft" type: self-excited oscillations 
emerge, with amplitudes that scale as y/k — kc, where 
kc is the critical coupling. These oscillations are of the 
combined qubit-oscillator system, and their time depen- 
dence can be characterized by a single phase variable. In 
addition, since the beam is the element that loses sta- 
bility, its dynamics enslaves those of the qubit, and the 
frequency of oscillations is equal to its renormalized fre- 
quency. In the terminology of nonlinear dynamics, the 
system is said to undergo a supercritical Hopf bifurca- 
tion. The quantitative details of this picture are provided 
in section (jlll C|) . 

A. Centered equations of motion 

The set of equations (|22l) is inhomogeneous, with a 
constant term in the right hand side. To proceed with 
the analysis, we first transform the coordinates to a set 
centered around an equilibrium point Sep. This equi- 
librium point is obtained by expansion in powers of fc, 
namely 

Sep = Sep,o + ^Sep^i -I- fc^SEP,2 + . . . , (28) 

where Sep,o is the single solution of the system when 
A; = 0. We find the following set of equations: 

Sc-{t) = -72Sc-(t) -i{s + k{a,{t) + «:(*))) s,_(t) 
+ iiA„s,;,(t) - iK{a^{t) + al{t)) (29a) 

Scz (t) = -JlScz (t) + I An {Sc- (t) - Sc+ (t)) (29b) 

dc(t) = -^-f,^acit) - iuj^adt) ~ -ikscz{t), (29c) 
where Sc{t) — S{t) — Sep, and 

k, = -CcqfcA„ {6 + 172) , (30) 

^ _ 7l'^z,eq 

We have also assumed that both the dissipation and the 
coupling are small, namely that 71 ~ 72 <C <^ — A„, 
7ni <^ Wni, and k <C cOm, and that uj^ is never significantly 
larger than 



6 



the Rabi frequency of the qubit. The sahent feature 
of (|29p is that now we see that the oscihator interacts 
with the qubit through a hnear driving term in ()29a|) . in 
addition to the nonhnear frequency-shifting term found 
in ((22a|) . 



B. Stability and renormalized mechanical 
eigenvalue 

Since the stabiUty of any dissipative dynamical system 
is determined solely by its linear part, we characterize the 
stability of (|29| by neglecting its nonlinearity, and find an 
expression for the renormalized mechanical eigenvalues, 
which we denote as X± — — ^7m i iuini- Linearization 
gives us an approximate equation of the form Sc = J'Sc, 
where 



J 





-72 — iS 


^ izA„ 









-72 + iS -\iAn 








-A„ -71 










-\ik 


"57m~«Wni 


V 





\ik 


-57m + ia;,„/ 



(31) 



is the Jacobian of the system. We are now facing the 
problem of finding the eigenvalues of J'. This can be 
achieved using perturbation theory, as shown in ap- 
pendix [X] We will follow, however, a different approach 
that yields identical results, yet has the benefit of pro- 
viding considerable physical insight and can be applied 
to general problems that involve linearly interacting sub- 
systems. 

We will focus our attention on the qubit and on the 
oscillator separately, regarding both of them as input- 
output systems. From the form of linear part of (I^Hl) . 
we can see that it is possible to regard adt) + a*(f) as 
the input of the qubit, and Scz(^) as the input of the 
osciUator. If we also regard Scz{t) as the output of the 
qubit and adt) + a*{t) as the output of the oscillator, 
we can decompose the matrix into two separate 3x3 and 
2x2 blocks. This simplifies the analysis considerably, and 
provides physical insight into the behavior of the qubit 
under linear driving of the type seen in (j29al) . As we 
will see, the stability of the system and the renormalized 
mechanical eigenvalues A± are determined mostly by the 
response of the qubit to the linear driving term in (I29ap — 
in particular, we will show in subsection IIII B 2l that 



fclm Xz, 



(32) 



where Xz is the qubit response function. We therefore 
begin by deriving an approximate expression for Xz and 
discussing its properties. 



1. Qubit response function 

To derive Xzi the response function of the qubit, we as- 
sume that the time dependence of the oscillator is a gen- 
eral complex sinusoid of the form ac{t) = aoe"*, where 



s is complex. Due to the linearity of the system, the 
steady state response of the z component of the qubit, 
Scz(i), will be given by Xz(s)aoe^*. This function, as well 
as the response functions of Sc- (t) and Sc+ (t) , are given 
by 




(3) 




(33) 



where J'^'^^ is the upper left 3x3 block of J', and is 
the 3x3 identity matrix. Invoking the assumptions on 
the values of the parameters made in section IIII Al and 
using = 72 — 71 , we have 

5CeqfcA^(s + 272) 



Xz{s) 



- 72 ■ 



■72 ■ 



(34) 



where 



71 = 71 + 



Al 

02 ^^P' 
"r 

A2 



72 = 72 ■ 



2rj| 



If, 



(35) 



472X + 7^A2 {A6' + Al) 



The location of the poles and zeros of (p4)) is correct to 
second order in the small parameters. A plot of the real 
and imaginary parts of Xzi—i'-^) can be found in Fig. [3l 
and a plot of Xz as a function of w and S is given in 
Fig. [2] Note that the correction for the Rabi frequency, 
fift — Or, is of second order in the small parameters, and 
thus can be neglected in what follows. Note that 71 and 
72 correspond to the approximate decay constants in the 
diagonal coordinate system of the qubit, provided that 
7i and 72 are small compared to JIr. This will be shown 
explicitly in appendix [Bl where the diagonal coordinates 
of the qubit will be found. 
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FIG. 2. (Color online) The response function Xzi^i^) as 
given in (|34|l . as a function of 5 and tJ. The vertical axis cor- 
responds to the absolute value, the color corresponds to the 
phase, and the contours correspond to the imaginary part. 
For all S < 0, the phase is positive, which implies that the 
qubit adds delay, thus decreasing the effective dissipation con- 
stant of the oscillator. (Recall that we are considering nega- 
tive frequencies.) For S > 0, the opposite is true. Parameters 
are given in Table HI 



8. Renormalization of mechanical eigenvalue 

To derive l|32|) . let us consider the qubit as a part of 
the full system. Since the Q-factor of the mechanical 
oscillator, Q = y^, is typically high (of the order ~ 10^), 
we can expect that any excitation at a frequency far from 
Win will not be sustained by the system, and that the 
decay or growth times of the oscillator are very long. 
This allows us to assume that the qubit is driven by a 
signal of the form Q;oe~*'^'"*-|-c.c., and due to linearity, the 
z-coniponent of the Bloch vector of the qubit responds 
as Scz{t) = Xz(~*'^m)aoe~*'^"* + c.c. Substituting this 
response in the equation for Uc in ()29c|) and keeping only 
resonant terms, we find that under this approximation 
we have 

ac{t) = -^7mac(i) - i^mOic[t)i (36) 
where 7„i and a)„i are given by 

3. Nearly adiabatic case 



(a) 




0.001 



0.015 r 
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(b) 

Re(X;(.-'<^)) 



0.01 0.02 0.05 0.10 0.20 0.50 1.00 2.00 



FIG. 3. (Color online) The imaginary (a) and real (b) parts of 
the qubit response function, Xz(— ii^), as given by (|34p in the 
blue-detuned {S < 0) regime, for different qubit decay times 
7i and 72. The blue dots in (a) corresponds to the maximum 
at oi = 71, and the red dots to the maximum ai u — Or. Solid 
line: 71 = 0.001, 72 = 0.01. Dashed line: 71 = 0.05. 72 = 0.1, 
Dot-dashed line: 71 = 0.1, 72 = 0.5. Other parameters are 
given in Tabled 



In this case we assume that Wm ^ JIr, and that Wm ~ 
71, 72, so that Wm is also a small parameter. Substituting 
s = — iwjn in (j34l) and keeping only small terms only up 
to second order, we find that 

(5CoqA:A2 (w„ 



2^72) 



'R V-^m -r 171) 

Substituting this into (P^ . we find that 

7m = 



and 



<5Ccqfc2A2 27271 



2^1 



(37) 



(38) 



(39) 



Note that because 71^ > 0, the sign of the correction to 7ni 
in ((38|) depends only on the sign of 5, and is negative for 
blue detuning and positive for red detuning, as expected. 

As we can see, the critical coupling kc^ that can be 
found by equating (p8|) to zero, depends strongly on 
the relation between the characteristic time scales of 
the qubit and the oscillator, and in particular on the 
ratio Wm/riR. The interaction is maximal at the reso- 
nance, Win/ilR — 1, and exhibits another smaller peak at 
Win = 7ij where 71 is defined in psp . Outside the range 
7i < Win < f^R, the interaction is weak. This dependence 
is illustrated in Fig. |3] which plots (l34t . From a control 
systems point of view, we can see that the qubit approx- 
imately responds as a harmonic oscillator with natural 
frequency JIr and decay constant 7in, connected in se- 
ries to a lag compensator with a much slower time scale. 
This implies that when (5 < 0, the qubit will reduce the 
effective dissipation coefficient for any set of values of the 
other parameters, and when (5 > 0, the opposite is true. 
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4- Resonance case 



1. Approximated equations of motion 



When — f^R, the qubit and the oscihator are res- 
onant, and the interaction is maximal. In this case, 
the correction of the mechanical dissipation coefficient, 
7ni — 7m, becomes of first order in the small parame- 
ters of the system, and in fact overwhelms 7ni for high 
quality oscillators. To derive an expression for Xz near 
resonance, we set 



The set of equations described above is given by 



Win = ^^R + CT, 

where ct is a small parameter. We then get 



(cr + i72) ' 



(40) 



(41) 



which gives 



and 



2r!R (7| + (t2) 



(42) 



(43) 



The general behavior of the correction to the mechan- 
ical dissipation coefficient, as a function of the DC and 
AC parts of the externally applied flux, eo^co f^nd So'/'ei, 
can be seen in Fig. |4l This result was obtained by su- 
perposing the corrections for different values of n, the 
multi-photon Rabi resonance, as defined in (|20|) . This 
result may be compared to the Landau-Zener interfer- 
ence diagrams, given in Refs. HH and 113. 



C. Hopf bifurcation at resonance regime 



Sc-{t) 
Scz{t) 

adt) 



: ~ (72 - ia) Sc- (t) 

ikKnac{t){scz{t) + Sz.oq), 
- -llScz{t) 



(44a) 
(44b) 

(44c) 



where 



An 

2nn 



a is the small detuning from resonance, defined in (|40p . 
and Sz^eq = <5f^RCoq- ThcsB equations are derived in ap- 
pendix |B1 using the procedure outlined above. Note that 
Sz.eq > for blue detuning — this corresponds to an effec- 
tive inversion of the energy levels of the qubit. 

Since these are envelope equations, we can obtain a 
detailed picture of the behavior of the system near res- 
onance by finding their equilibrium points. Specifically, 
we will show that for 7ni > 0, where 7ni is as defined 
in (I42p . the system will have a single stable equilibrium 
at the origin, while for 7m < 0, the origin will lose its sta- 
bility and a new set of equilibrium points will emerge. In 
the approximated system (j44p . this change is a supercrit- 
ical pitchfork bifurcation, while for the original equations 
in (|29p , it corresponds to a supercritical Hopf bifurcation. 



2. Limit cycle solution 

To find the non-trivial equilibrium points of (|44l) . we 
begin by assuming that during the limit cycle, the phase 
difference between Oc and Scz remains constant, which is 
required by the periodicity of the limit cycle. Therefore, 
we can set: 



In this section, we no longer neglect the nonlinear part 
of ([29|l . We also continue assuming that <5 < so that the 
system is blue-detuned. We will derive approximate an- 
alytical expressions for the limit cycle amplitude and fre- 
quency, under the assumption that the system is approx- 
imately resonant, i.e. that Wm — ^k, using the method of 
averaging (chapter 11 of Ref. IsgI ) . This is accomplished 
by first diagonalizing Mq, the large part of M in (pS)) de- 
fined in (|B1[) , and then transforming to a frame rotating 
at Wm- This will essentially give us an equation for the 
envelopes of components of the qubit on the Bloch sphere 
and of the amplitude of the oscillator, which are slowly 
varying. The equilibrium points of this set of equations, 
then, will give the amplitude of the limit cycle, while the 
equation for the phase of these envelopes will provide the 
frequency of oscillations. 



ttc = Ta exp [iuSat) 



(45) 



Ts exp 



where r^, r^, Wq and /(cr) are real, and the form for 
the phase of Sc- was chosen so that /(cr) = for cr = 0. 
When the system is in a limit cycle state, both the radii 
and Ta and Scz are constant in time, and we have assumed 
that uJa is constant as well. Therefore, substitution of 
these definitions into gives a set of equations for the 
equilibrium point that can be solved exactly. This set 
has a nontrivial solution only for k > kc, where 



kn — 



'27mf^R (72^+^2) 



'5Coq72A2 



(46) 
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(a) (b) 




20 40 60 80 100 20 40 60 80 100 



FIG. 4. (Color online) Correction of dissipation coefficient as a function of external flux driving amplitude, eo(t>ci/h, and 
detuning, eo(l)co/h. The blue areas corresponds to a decrease, and the red to an increase, in the effective mechanical dissipation 
coefficient. The color intensities are scaled by the extremal values of the correction. Panel (a) corresponds to the resonance 
regime, with iUm — 1.28A, and panel (b) corresponds to the nearly adiabatic regime, with = 0.14A. The values of the 
correction are given in arbitrary units. Parameters used are: 71 — 0.014ajd, 72 ~ 0.714ajd, k = O.OOlSoJd, o-z,cq ~ —1, A — O.lcJd, 
and Q — 10^. This plot may be compared to Refs. [31I and [34. 



and we have neglected 7m in comparison with 72. Note 
that this is the same critical coupling as that given 
by (H^ . (The assumption that 72 ^ 7m was made in sec- 
tion |IITB] when we assumed that the input to the qubit is 
a non-damped sinusoid.) The nontrivial solution is given 
by: 




(47) 



7m cr 



272 + 7n 



/(ct) = arctan 



2a 



272 + 7n 



This solution provides a complete characterization of the 
limit cycle near resonance. We can see that the system 
has a limit cycle for k > kc, with amplitudes proportional 
to \/k — kc, as is the case with non-degenerate supercrit- 
ical Hopf bifurcations. Furthermore, the frequency of 
the limit cycle solution undergoes a shift proportional 
to a, as can be seen in the equation for Ua in (1471) . In 
Fig. [5] we can see the bifurcation curves for the system 
near resonance, which were calculated using numerical 
continuatioii^. 



In the nearly adiabatic regime the interaction is of 
second order in k. It turns out that bifurcation anal- 
ysis of the approximated equations of the form pi)) for 
this regime is analytically intractable. Numerical analy- 
sis, however, indicates that the Hopf bifurcation in this 
regime is supercritical as well, with behavior similar to 
that seen in the resonance regime. 



IV. DISCUSSION AND CONCLUSIONS 

We have characterized the stability properties of a flux 
qubit coupled to a mechanical oscillator and its nonlinear 
behavior near resonance. We have found that the interac- 
tion between the qubit and the oscillator is typically non- 



Parameter 



Value 



71 
72 

S 
A 

k 



O.OIS^R 

-0.707J7r 
0.707!^R 
0.013f7R 
-1 
10^ 



TABLE I. Numerical parameters used in the generation of 
Figures [1 [2] and El 
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0.02 



0.04 0.06 



0.02 



0.04 0.06 



0.08 




0.08 



0.06 
0.04 
0.02 
0.00 
-0.02 
-0.04 



0.02 



0.04 0.06 



0.02 



0.04 0.06 



0.08 




0.08 



FIG. 5. (Color online) Bifurcation curves for the system near resonance, with ujm = LIHr. The blue and dashed red lines 
corresponds to stable and non-stable equilibrium points, respectively. The green curve corresponds to a stable limit cycle, 
where the top curves corresponds to the maximal value of the variable during the limit cycle, and the bottom curve corresponds 
to the minimal value. Numerical values used in the simulation are given in table HI 



negligible only over a range of oscillator frequencies that 
begins with Wm — 71, which we denoted as the nearly- 
adiabatic regime, and ends with oj^ ~ JIr, which we 
denoted as the resonance regime. The interaction near 
the resonance regime is significantly stronger. We have 
found that sign of the correction factor for the mechani- 
cal dissipation coefficient 7ni depends only on sign of the 
detuning parameter S. 

We have also found that near resonance, when the crit- 
ical coupling is exceeded, the system undergoes a super- 
critical Hopf bifurcation. Numerical analysis indicates 
that the bifurcation at the nearly adiabatic regime is su- 
percritical as well. 

The method used in the stability analysis of this sys- 
tem, namely the decomposition into interacting subsys- 
tems, is general in nature and can be applied to any set 
of weakly interacting subsystems, yielding considerable 
physical insight into their behavior. 

The limit cycle behavior of this system suggests a pos- 
sible scheme for the preparation of non-classical entan- 
gled states of the qubit and oscillator. Since the limit 
cycle dynamics can be described by a single phase vari- 
able, if the system is cooled to its ground state so that 
thermal noise is negligible and then rapidly brought to a 



limit cycle state, this phase variable can be expected to 
be found in a superposition state. A more quantitative 
analysis of this point is left to subsequent articles. 
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Appendix A: Perturbation analysis of mechanical 
eigenvalue 

The main results regarding the renormalized mechani- 
cal eigenvalues A±, jMl), (EH), dMl) and (gS]), can also be 
obtained by using perturbation analysis for J^, under the 
assumption that fc is a small parameter. Let us denote 

J = Jo + kJi, (Al) 
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where kj'i is the part of J7 proportional to k. We also de- 
note Ao+ — — ^7m + i'^m as the mechanical eigenvalue of 
Jo, and as its corresponding eigenvector. We now de- 
note 7?. as the inverse of the upper 4x4 part of A+ /s — Jo, 
namely the matrix satisfying 

(X+h -Jo)n^n (A+/5 - Jo) - /5 - v„,vT , (A2) 

where I5 is the 5x5 identity matrix. It is then possible 
to show that A-|_ can be expanded in a power series in fc, 
where 



Ao+ + k^^lJk^m + k^^rlJlJk'Rw^ + Oi^). (A3) 

The term proportional to k is equal to zero, and the term 
proportional to k^ gives the approximate results obtained 
above. 



Appendix B: Approximated Equations for resonance 
regime 

To derive the approximated equations for the en- 
velopes (|44p . we begin by diagonalizing the large part 
of as explained above. Let us define M — Mq + Mi. 
In the basis defined by ((27|) . 



and 



-8 


-A„ 


iA„ 
5 -iA„ 
A„ 


^ 


(Bl) 





-CJin 
/ 





Mo = i 



is the large part of M, and therefore of (|26|) . and Mi 
contains the smaller dissipative terms on the diagonal. 
The diagonalizing matrix of the upper 3x3 block of this 
Mo, which we denote as Mq , is given by 



)(3) 







1^ 


\ 




:! 




I Or + 




IOr 


\ 


1 ^ - 




\ 


1 ^ + 






I Or 




I, Or ^ 




A„ 






A„ 






Or 






Or 





20n 

A„ 
20n 
b 



(B2) 



-il^R 0^ 

il^R 
0, 



Next, we note that close to resonance, a 
small parameter. Therefore, 



(B3) 



f^R is a 



i^(3) 















can be separated to a large part -Do, proportional to 
Win, and to a small part proportional to a. We can 
now define a set of envelope coordinates Sc(t) such that 
Sc(t) = Qe^«*Sc(t). Explicitly: 



Sc(i) - (s,_(i),Se+(t),Sc.(t),ac(i),a:(t))^ 



(B4) 



and 





' 








where 1-2 is the 2x2 identity matrix. When substitut- 
ing (|B4[) into we obtain the equation 



Sc(t) = e-^"*Q-i (I?, + Ml) Qe^«*Sc(t) 
+ /ce-^«*Q-ig(Qe^«*S,(i)) , 



(B5) 



Since all terms on the right hand side of this equation 
are small and e^°* contains a single frequency, we 
can now perform the averaging approximation. This is 
achieved by integrating (jB5p over a single period and 
keeping Sc(i) constant: 



1 



2T:U)m^ Jo 
k 



27ra;in^ Jo 



e-^o* Q-ig I^QgA,* s,{t)j dt'. (B6) 



The integration will remove all terms in (|B5P propor- 
tional to e'""* and its harmonics, and is in fact equiv- 
alent to performing an RWA. Equation (IB6p is equal to 
the set (Ii4l). 
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